lambda_seq <- 100*seq(1,.505,-.005)^14
i = 0
lasso.results <- list()
for (year in start.year:end.year) {
  i = i + 1
  print(year)
  # Fit LASSO
  lasso.results[[i]] <- list()
  dta.past <- dta[dta[,t.var]< year,]
  N <- length(dta.past[,dv])
  cv_results <- foreach (cv=1:cv.runs, 
                         .combine='rbind',
                         .packages = c("glmnet","pROC")) %dopar% {   # Loop over CV runs
                           set.seed(52*cv) 
                           shuffle <- sample(1:N,
                                             N,
                                             replace=FALSE)
                           prediction <- matrix(NA,nrow=N,ncol=length(lambda_seq))
                           sampleSplit = c(0,round((N/5*c(1:5))))
                           for (f in 1:5) {
                             # Fit lasso
                             cv.fit <- glmnet(x = as.matrix(dta.past[-shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                                          rhs]),
                                              y = as.matrix(dta.past[-shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                                          dv]),
                                              alpha = lasso_count_params$alpha,
                                              family = lasso_count_params$family,
                                              lambda = lambda_seq)
                             prediction[shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                       1:length(cv.fit$lambda)] <- predict(cv.fit,
                                                                           newx = as.matrix(dta.past[shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                                                                          rhs]),
                                                                           type = "response")
                           }
                           cv_results <- c("CV Run"=0,
                                           "lambda"=0,
                                           "R2"=0)
                           # Get error rates for each lambda used
                           for (l in seq(10,90,10)) {
                             if (sum(is.na(prediction[,l]))==0) {
                               R2 <- Rdev(dta.past[,dv],
                                          prediction[,l])
                               if (R2>cv_results["R2"]) {
                                 cv_results <- c("CV Run"=cv,
                                                 "lambda"=lambda_seq[l],
                                                 "R2" = R2)
                                 best.l <- l
                               }
                             }
                           }
                           for (l in seq(best.l-9,best.l+10,1)) {
                             if (sum(is.na(prediction[,l]))==0) {
                               R2 <- Rdev(dta.past[,dv],
                                          prediction[,l])
                               if (R2>cv_results["R2"]) {
                                 cv_results <- c("CV Run"=cv,
                                                 "lambda"=lambda_seq[l],
                                                 "R2" = R2)
                                 best.l <- l
                               }
                             }
                           }
                           return(cv_results)
                         } 
  hours <- floor((proc.time()[3]-start.time)/3600)
  mins <- floor((proc.time()[3]-start.time)/60) - hours*60
  secs <- floor((proc.time()[3]-start.time)) - hours*3600 - mins*60
  
  print(paste(hours,
              "h",
              mins,
              "m",
              secs,
              "s elapsed",
              sep = ""))
  
  optim.lambda <- mean(cv_results[,"lambda"]) 
  lasso.results[[i]]$mod <- glmnet(x = as.matrix(dta[dta[,t.var]< year,
                                                          rhs]),
                                   y = as.matrix(dta[dta[,t.var]< year,
                                                          dv]),
                                   alpha = lasso_count_params$alpha,
                                   family = lasso_count_params$family,
                                   lambda = seq(991,1,-10)*optim.lambda
  )
  opt.lam <- 100
  while (is.na(lasso.results[[i]]$mod$lambda[opt.lam])) {
    opt.lam = opt.lam - 1
  }
  lasso.results[[i]]$lambda <- optim.lambda
  lasso.results[[i]]$fit.oos <- predict(lasso.results[[i]]$mod,
                                        newx = as.matrix(dta[dta[,t.var] == year,
                                                                  rhs]),
                                        type = "response")[,opt.lam]
  lasso.results[[i]]$actual.oos <- dta[dta[,t.var] == year,
                                            dv]
}
save(lasso.results,
     file=paste(modeldir,"/",
                country,
                "_lasso_",
                "count",
                "_",
                rhs.group,fileext,
                ".RData",
                sep = ""))